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Abstract 

Genetic switch systems with mutual repression of two transcription factors are studied using 
deterministic methods (rate equations) and stochastic methods (the master equation and Monte 
Carlo simulations). These systems exhibit bistability, namely two stable states such that spon- 
taneous transitions between them are rare. Induced transitions may take place as a result of an 
external stimulus. We study several variants of the genetic switch and examine the effects of 
cooperative binding, exclusive binding, protein-protein interactions and degradation of bound re- 
pressors. We identify the range of parameters in which bistability takes place, enabling the system 
to function as a switch. Numerous studies have concluded that cooperative binding is a necessary 
condition for the emergence of bistability in these systems. We show that a suitable combination 
of network structure and stochastic effects gives rise to bistability even without cooperative bind- 
ing. The average time between spontaneous transitions is evaluated as a function of the biological 
parameters. 

PACS numbers: 87.10. +e,87.16.-b 
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I. INTRODUCTION 



Recent advances in quantitative measurements of gene expression at the single-cell level 
1,0] have brought new insight on the importance of stochastic fluctuations in genetic circuits 
3|. The role of fluctuations is enhanced due to the discrete nature of the transcription 
factors and their binding sites, which may appear in low copy numbers . As a result, 
populations of genetically identical cells may show significant variability. Stochastic behavior 
may invoke oscillations [Q] and spatio-temporal patterns which are unaccounted for 
by macroscopic chemical rate equations. Genetic circuits with feedback mechanisms may 
exhibit bistability, namely, two distinct stable states which can be switched by an external 
signal Q]. A low rate of spontaneous switching events may also take place. To qualify 
as a switch, this rate must be much lower than the rates of the relevant processes in the 
cell, namely transcription, translation, and degradation of transcription factors. Genetic 
switches, such as the phage A switch, give rise to different cell fates joj]. In this switch, A 
phages infect E. coli bacteria and can exist in two exclusive states, one called lysogeny and 
the other called lysis. When the phage enters its host, it integrates itself into the host's 
DNA and is duplicated by cell division. It codes for proteins that can identify stress in the 
host cell. In case of stress, the phage transforms into the lysis state. In this state, it kills 
the host cell, using its DNA to produce many copies of the phage, which are released and 
later infect other cells. Other switch circuits exist in the metabolic systems of cells. These 

n 

switches determine which type of sugar the cell will digest . The genetic switch may also 
serve as a memory unit of the cell, and help determine its fate during cell differentiation. 

Recent advances enable the construction of genetic circuits with desired properties, that 
are determined by the network architecture. These networks are constructed from available 
components, namely genes and promoters. They do not require the manipulation of the 
structure of proteins and other regulatory elements at the molecular level. These genes and 
promoters are often inserted into plasmids rather than on the chromosome. A synthetic 
toggle switch, that consists of two repressible promoters with mutual ne gati ve regulation, 
was constructed in E. coli and the conditions for bistability were examined ll(| . The switch- 
ing between its two states was demonstrated using chemical and thermal induction. More 
recently, such circuit was found to exist in a natural system in which two mutual repressors 
regulate the differentiation of myeloid progenitors into either macrophages or neutrophils 
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In this paper we analyze the genetic toggle switch using deterministic and stochastic 
methods. In this simple genetic circuit, two proteins, A and B, negatively regulate each 
other's synthesis. The regulation is performed at the transcription level, namely the pro- 
duction of protein A is negatively regulated by protein B, through binding of n copies of B 
to the A promoter (and vice versa). This process can be modeled by a Hill function, which 
reduces the production rate of A by a factor of 1 + where [B]_is the concentration of 

B proteins in the cell, is a parameter and n is the Hill coefficient 13|]. In case that n = 1, 
the binding of a single protein is sufficient in order to perform the negative regulation, while 
for n > 1 the cooperative binding of two or more proteins is required. 

One may expect this circuit to function as a switch, with two stable states, one dominated 
by A proteins and the other dominated by B proteins. When the population of A proteins 
is larger than the population of B proteins, the A proteins suppress the production of B 
proteins. Under these conditions, the production of A proteins will not be suppressed much 
by the small B population. Therefore, the system approaches a state rich in A proteins and 
poor in B proteins. Similarly, the system may approach a state rich in B proteins and poor 
in A proteins. 

To qualify as a switch, the system should be bistable. In the deterministic description, 
bistability is defined as the existence of two stable steady state solutions of the rate equations. 
This description does not account for the possibility of spontaneous transitions between the 
two states. In the stochastic description, spontaneous transitions do take place. Therefore, 
the condition for bistability is that the rate of spontaneous switching events (due to random 
fluctuations rather than an external signal) is much lower than the rates of all other relevant 
processes in the system. 

Rate equations provide the average concentrations of A and B proteins in a population 
of cells. In these equations, bistability emerges at a bifurcation point, where two stable 
states emerge. Rate equations do not include fluctuations and do not account for the pos- 
sibility of spontaneous transitions between the two states. The master equation provides 
the probability distribution of the populations of A and B proteins. The two bistable states 
appear as two distinct peaks in this distribution. Monte Carlo simulations enable to follow 
the fluctuations in a single cell and to evaluate the rate of spontaneous switching events. 

We examine the conditions for the system to become a switch, and calculate the rate 
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of spontaneous transitions between its two states. This is done for several variants of the 
toggle switch. In particular, we focus on switch systems in which the repression in done 
without cooperative binding (namely, n = 1). Numerous studies have concluded, using rate 
p.^^th^ ^perat.ve bmdmg .s a necessary cond.t.on for the emergence of HstaHhty 



[111, 



14,^ 
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17(1 . Below we show, using a combination of deterministic and stochastic 
simulation methods, that this is not the case, namely a bistable switch can exist even in the 
absence of cooperative binding. In particular, we show that bound-repressor degradation 
(BRD) and protein-protein interactions (PPI) give rise to bistability, without cooperative 
binding, even at the level of rate equations. These results are confirmed by stochastic 
simulations using the master equation and Monte Carlo methods. We also consider the 
exclusive switch, in which the A and B repressors cannot be bound simultaneously due to 
overlap between their promoter sites. This system exhibits bistability only when stochastic 
fluctuations are taken into account. The rate of spontaneous transitions between the two 
states is calculated as a function of the biological parameters. 

The paper is organized as follows. In Sec. [Ill we consider the basic version called the 
general switch. Several variants of this circuit are considered in the Sections that follow. 
The exclusive switch is studied in Sec. Illlt the BRD switch is considered in Sec. IIVI and 
the PPI switch is analyzed in Sec. |Vl The effects of cooperative binding are studied in Sec. 
IVII The response of toggle switch systems to external signals is examined in Sec. I VI II The 
results are discussed in Sec. IVIIII and summarized in Sec. HXl 

II. THE GENERAL SWITCH (WITHOUT COOPERATIVE BINDING) 

The general switch consists of two transcription factors, A and B, that negatively regulate 



each other's synthesis I3, ll5|. A schematic description of this circuit is given in Fig. [T](a). 
The regulation is done by the binding of a protein to the promoter site of the other gene, 
blocking the access of the RNA polymerase and suppressing the transcription process. In 
this circuit there is no cooperative binding, namely the regulation process is performed by 
a single bound protein. 

The concentrations of free A and B proteins in the cell are denoted by [A] and [B], 
respectively (by concentration we mean the average copy number of proteins per cell). The 
copy numbers of the bound proteins, are denoted by [ta] and [tb], where va is a bound A 
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protein that monitors the production of B, while is a bound B protein that monitors the 
production of A. Note that there is at most one bound repressor of each type at any given 
time, and thus < ta^tb < 1- For simphcity, we ignore the mRNA level and combine the 
processes of transcription and translation as a single step of synthesis [l8|. 

The maximal production rate of protein X is denoted by gx (s~^), X = A,B. The 
degradation rate of protein X is given by dx (s~^). While the structure of the circuits 
studied here is symmetric, the rate constants can be different for A and B. However, for 
simplicity we use symmetric parameters, i.e. g = Qa = 9b and d = = ds- The binding 
rate of proteins to the promoter is denoted by ao (s~^) and the dissociation rate by ai (s~^). 

A. Rate Equations 

The dynamics of the general switch circuit is described by the rate equations [lol, Q 



[A] = gAil - [tb]) - dA[A] - ao[A] (1 - M) + aiM 

[B] = gB{l - M) - dB[B] - ao[B] (1 - [r^]) + a,[rB] 
[ta] = ao[A]{l- M) -ttiM 

[tb] = ao[B]{l-[rB])-a,[rB]. (1) 

It is commonly assumed that the binding-unbinding processes are much faster than other 
processes in the circuit, namely ao,ai ^ dx,gx- This means that the relaxation times of 
[rx] are much shorter than other relaxation times in the circuit. Under this assumption, one 
can take the time derivatives of [rx] to zero, even if the system is away from steady state. 
This brings the rate equations to the standard Michaelis-Menten form 



where symmetric parameters are used, and = ao/cn is the repression strength. For a given 
population of free X repressors, the parameter k controls the value of [rx]- The limit of 
weak repression, [rx] ^ 1, is obtained when k[X] ^ 1, while the limit of strong repression, 
[rx] — 1, is obtained for k[X] ^ 1. 
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The meaning of bistability at the level of rate equations is that at steady state the 
equations exhibit two distinct positive solutions. In this particular class of circuits, one 
solution is dominated by A proteins and the other is dominated by B proteins. Starting 
from any initial state, the system will converge to one of these solutions. The solutions are 
stable, so the possibility of spontaneous transitions, induced by stochastic fluctuations, is 
not included in the rate equation description. 

The steady state solutions of Eqs. ([1]) and ([2]) are identical. We will now show that 
these equations have only one positive steady-state solution. To this end, we first take 
[A] = [B] = in Eq. ([2]). We multiply each equation by the denominator of the Hill 
function that appears in it. We obtain: 



g - d[A] - kd[A][B] = 

g-d[B] - kd[A][B] = 0. (3) 

Subtracting one equation from the other we get d{[A] — [B]) = and therefore [A] must be 
equal to [B] at steady state. The steady state values of [A] and [B] can be easily found. 
Inserting [A] = [B] into Eq. ([3]) we obtain a quadratic equation whose only positive solution 
is 

-1 + Jl + Akg/d 
[A] = [B] = \^ (4) 

Standard linear stability analysis shows that this solution is always stable. 

As a result, we conclude that at the level of rate equations the general switch, with- 
out cooperative binding, does not exhibit bistability. In Sec. I VI I we consider the case of 
cooperative binding, where the rate equations do exhibit bistability. 



B. Master Equation 



In order to account for stochastic effects and to obtain insight on the reason that this 
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24| . In this case. 



system is not bistable, the master equation approach is applied [3|, l 
we consider the probability distribution function P{Na, Nb, ta, tb). It is the probability for 
a cell to include Nx copies of free protein X and rx copies of the bound X repressor, where 
Nx = 0, 1, 2, . . ., and rx = 0, 1. The master equation for the general switch takes the form 



P{Na, Nb, rA, rn) = gx^r^AP^NA - 1, A^^, r^, tb) - P{Na, Nb, ta, tb)] 

+gBSr^,o[P{NA, Nb - 1, rA, Tb) - P{Na, Nb, ta, Tb)] 

+dA[{NA + 1)P{Na + 1, Nb, t-a, r^) - NaP{Na, Nb, ta, tb)] 

+dj,[{NB + 1)P{Na, Nb + 1, Ta, tb) - NbP{Na, Nb, ta, tb)] 

+ao[(A^A + V)br^,xP{NA + 1, Nb, ta - 1, tb) - NASr^^P^NA, Nb, ta, tb)] 

+aQ[{NB + l)5r^,iP{NA, Nb + 1, ta, - 1) - iVB5,^,oP(iVA, Nb, ta, tb)] 

+ai[5r^^QP{NA - 1, Nb, ta + 1, r^) - 5r^,iP(A^A, Nb, ta, tb)] 

+ai [5r^flP{NA, Nb - 1, ta, tb + I)- 5r^,iP{NA, Nb, ta, tb)], (5) 

where 5i^j = 1 for i = j and otherwise. The gx terms account for the production of proteins. 
The dx terms account for the degradation of free proteins, while the ao (ai) terms describe 
the binding (unbinding) of proteins to (from) the promoter site. In numerical integration, 
the master equation must be truncated in order to keep the number of equations finite. This 
is done by setting suitable upper cutoffs, N^^^ and A^^'^^, on the populations sizes of free 
proteins. In order to maintain the accuracy of the calculations, the probability of population 
sizes beyond the cutoffs must be sufficiently small. 

The master equation has a single steady state solution, which is always stable 251]. The 
criterion for bistability is that the steady state solution P{Na, NB,rA,rB) exhibits two 
distinct regions (peaks) of high probabilities, separated by a gap in which the probabilities 
are very small. These two regions correspond to the two states in which the system is likely 
to be. If the transition rate between the peaks is small enough, the system is indeed a 
bistable switch. Note, that in this case, averages of the form 

{Nx)= E E E T.^^PiNA,NB,rA,rB), (6) 

Na=0 Nb=0 r'A=0 rB=0 

where X = A, B, do not reflect the complex structure of the probability distribution. These 
can be considered as averages over many cells, some dominated by A and others dominated 
by B proteins, such that the total populations of the two species are about the same. 

To examine the existence of bistability we consider the marginal probability distribution 
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P{Na,Nb) = E E P{NA,NB,rA,rB). (7) 

'■a=0 r-B=0 

This probability distribution was calculated for a broad range of parameters. Two represen- 
tative examples are shown in Fig. [2l 

Under conditions of weak repression (small k), P{Na, Nb) exhibits a single peak for which 
A^^ ^ Nb ^ g/d [Fig. [21(a)], in agreement with the rate equations. This is due to the fact 
that the repression is weak, and the A and B populations are almost uncorrelated. In this 
case, the cell will contain roughly the same amount of A and B proteins. 

For strong repression, the distribution P{Na, Nb) exhibits a peak dominated by A pro- 
teins and a peak dominated by B proteins, as expected for a bistable system. However, a 
third peak appears near the origin, in which both populations of free proteins are suppressed 
[Fig. El^b)]. This peak represents a dead-lock situation, caused by the fact that both A and B 
repressors can be bound simultaneously, each bringing to a halt the production of the other 
specie. The third peak provides a corridor through which the probability can flow between 
the other two peaks. As a result, the system can quickly switch between the A-dominated 
and the S-dominated states. 

In addition to the solution of the master equation, Monte Carlo simulations have been 
performed. In these simulations one can follow the time evolution of the populations of 
free and bound proteins in a single cell. In Fig. [3]^a) we present the population sizes of 
free proteins vs. time for the general switch. It is clear that the cell can indeed be in one 
of three states: a state rich in A, a state rich in B and a state in which both proteins 
are in very low copy numbers. We conclude that a necessary condition for the system to 
become a switch is to prevent this dead-lock situation in which both protein populations 
are suppressed simultaneously. Below we present several variants of the circuit in which the 
third peak is suppressed, giving rise to a bistable switch. 



III. THE EXCLUSIVE SWITCH 

The first variant we consider is the exclusive switch, depicted in Fig. [T](b). In this circuit 
there is an overlap between the promoters of A and B. As a result, there is no room for both 
A and B proteins to be bound simultaneously. Exclusive binding is encountered in nature, 
for example, in the lysis- lysogeny switch of phage A 
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It was shown that in presence of cooperative binding, the exclusive switch is more stable 
than the general switch Q, Q- This is because in the exclusive switch the access of 
the minority proteins to the promoter site is blocked by the dominant proteins. Here we 
show that in the exclusive switch, stochastic effects give rise to bistability even without 
cooperativity between the transcription factors. The dead-lock situation in prevented in 
this case, since A and B repressors cannot be bound simultaneously. 



A. Rate Equations 

To model the exclusive switch, recall that the variable [ta] {[i^b]) is actually the fraction 
of time in which the promoter is occupied by a bound A (B) protein [l^. The fraction of 
time in which the promoter is vacant is thus 1 — [r^] — [r^]. Incorporating this into Eq. ([1]) 
gives rise to the following modification: in the ao terms, each appearance of [ta] or [tb] 
should be replaced by [r^] + [r^]. With this modification, the rate equations take the form 



= gil - - d[A] - ao[A] (1 - M - [tb]) + «iM 
[B] = g{l - M) - d[B] - ao[B] (1 - [ta] - [tb]) + «iM 
[ta] = ao[A]{l- [ta] - [tb]) - "iM 

[rs] = ao[B] (1 - M - [r^]) - a,[rB]. (8) 

Under steady state conditions, the rate equations can be reduced to the Michaelis-Menten 
form 



where, as before, k = a^jax. We will now show that even for the case of the exclusive switch, 
the rate equations still exhibit a single solution, thus there is no bistability. This is done by 
taking \A\ = [5] = and getting rid of the denominators, by repeated multiplications. The 
resulting equations are 



g^{kg-d)\A\~kd\A\{\A\-r\B\) = 
9 



g + ^kg-d)[B]-kd[B]{[A] + [B]) 
By subtraction of one equation from the other, we find that 



= 0. 



(10) 



{kg~d-kd{[A] + [Bm[A]-[B])=Q. 



(11) 



The positive, symmetric solution, [A\ = [B], is given by 



[A] 



{kg -d) + ^{kg + dy + Akgd 



(12) 



Akd 



The other, non-symmetric solution, given by 



kg-d-kd{\A] + [B]) = 



(13) 



is inconsistent with Eq. flTUl) unless g = 0, namely there is no production of A and B proteins, 
which immediately leads to [A] = [B] = 0. Under these conditions the solution of Eq. f|T3|) 
is [A] + [B] = —l/k, which requires a negative population size and thus makes no physical 
sense. Therefore, the only solution for g > is the symmetric solution, [A] = [B]. Thus, 
the rate equations do not support a bistable solution for the exclusive switch for any choice 
of the parameters. 

B. Master Equation 

To account for the effects of fluctuations, we now describe the exclusive switch using 
the master equation. It is similar to to master equation for the general switch given by 
Eq. ([5]), except for the following modifications: (a) In the ao and ai terms, each time 
^rA,j (^rsj), j = 0,1, appears it should be multiplied by ^r^.O; (^r4,o); (b) The constraint 
P{Na, Nb, 1, 1) = should be imposed. Implementing these changes we obtain the following 
equation: 



P{Na, Nb, rA, r^) = gASr^APi^A - 1, Nb, ta, rs) - P{Na, Nb, ta, tb)] 
+gBSrA,o[P{NA, Nb - 1, r^, r^) - P{Na, Nb, ta, tb)] 
+dA[{NA + l)P(iVA + 1, Nb, ta, tb) - NaP{Na, Nb, ta, tb)] 
+ciB[(iVB + ^)P{Na, Nb + 1, TA, Tb) - NbP{Na, Nb, rA, tb)] 
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+aoSr^fl[{NA + l)Sr^,iP{NA + 1, A^^, rA - 1, r^) - NASr^flP{NA, Nb, rA, rs)] 
+ao5,^,o[(A^B + l)5r^,iP{NA, Nb + 1, r^, - 1) - NB5r^,oPiNA, Nb, ta, tb)] 
+ai [5.^,oP(A^^ - 1, Nb, ta + 1, r^) - 5.^,iP(Ar^, A^^, ta, tb)] 

+ai [^.^^oi'lA^A, - 1, TA, TB + l)- 5r^,iP{NA, Nb, ta, Tb)]. (14) 

For the exclusive switch, as for the general switch, under conditions of weak repression, 
P{Na,Nb) exhibits a single peak [Fig. Hl^a)] that satisfies Na ~ Nb ^ g/d. However, as 
the repression strength increases two distinct peaks begin to form. For intermediate values 
of k these peaks are still connected, by a corridor of non- vanishing probabilities [Fig. Hl^b)]. 
Monte Carlo simulations show that for intermediate values of k, the system indeed exhibits 
two states, one rich in A and the other rich in B, but rapid transitions occur between them. 

For strong repression, the distribution P{Na, Nb) exhibits two peaks which are separated 
by a region with vanishing probabilities [Fig. |4](c)]. In one peak the A population is sup- 
pressed, while in the other peak the B population is suppressed, as expected for a bistable 
system. The average population of the dominant protein specie in each peak is (A^x) ^ g/d, 
while the population of the suppressed specie is (A^x) ~ 0. Monte Carlo simulations show 
that in this case the average time between spontaneous transitions is much longer. The 
typical switching time for the case shown in Fig. [3]^b) is around 10^ seconds. It is much 
longer than the time-scales of the transcription, translation and degradation processes. It is 
also longer than the time between cell divisions which is of the order of 10^ — 10^ seconds. 
The Monte Carlo results clearly show a large number of failed attempts in which a protein 
of the minority specie binds to the promoter and then unbinds again, without causing the 
system to flip. 

C. Analysis of Switching Times 

To evaluate the switching times we performed the following procedure. We initialized the 
master equation in a state which is completely dominated by A proteins, namely, P{Na = 
[g/d\,NB = 0, = 0, = 0) = 1 (where [ J represents the integer part), and all other 
probabilities vanish. The master equation was then integrated numerically and P{Na, Nb) 
was calculated as a function of time. The function f(t) = P{Na > Nb) — P{Na < Nb) 
was found to decay exponentially from its initial value, /(O) = 1, to zero, according to 
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/(t) = exp(— t/r). The time constant r is defined as the switching time 26|. Its inverse, 
r~^, is referred to as the switching rate. 

Using this procedure, we examined the dependence of the switching time, r on the protein 
synthesis rate, g [Fig. El^a)], the degradation rate, d [Fig. E](b)], and the repression strength, 
k [Fig. E](c)]. In the parameter range in which bistability takes place, we obtain that (a) 
T g, (b) r ~ l/d"^ and (c) t k. Concernin g F ig. [5]^c), note that system exhibits 



bistability only in the regime in which k is large [27||. For /c<l,r~100 — 1000 (s), 
which is the typical time-scale of other processes in the cell. Only for A; > 10, r becomes 
significantly larger than the time scales of other processes, and the system can function as 
a stable switch. The scaling properties of the switching time can be summarized by 

T — . (15) 

This result can be reproduced by a simple argument. Consider an initial state in which the 
system is dominated by A proteins, while the population of B proteins is suppressed, namely 
[A] ^ [B]. In this situation the promoter site is occupied by an A protein during most of 
the time. In order that the switch will flip, the bound A protein must unbind (at rate ai). 
Then, a B protein (rather than an A protein) should bind to the promoter. The probability 
for this to happen is ~ [5] /[A]. This B protein should remain bound long enough in order 
to build up a sufficiently large population of B proteins. On average, the B protein stays 
bound l/tti (s), during which g/ai proteins of type B are produced. After the B repressor 
will unbind, the probability that the next protein that binds will be of type B rather than A, 
is thus ~ {g/ai)/[A] (neglecting the degradation of A proteins, because ai d). Following 
this argument, the switching rate is given by 



[B] [B] 
[A] aM] ^[A] 
From the Michaelis-Menten equations we obtain that 

M = ^^.J- (17) 

[A] l + k[A] k[Ay ^ ' 

since for strong repression k\A\ ^ 1. Inserting this result into Eq. (fT6l) and using and 
\A\ g/dwe find that 
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(18) 

This result can be considered as the leading term in the expansion of r in powers of g, d and 
k. This leading term turns out to provide a very good approximation to simulation results. 
For example for g = 0.2, d = 0.005, ao = 0.2, and ai = 0.01 (s"^) we get r = 1.6 ■ 10^ (s), 
which agrees perfectly with the results of Monte Carlo simulations. 

From Eq. (fT6!) . and from the fact that the average copy number of the dominant specie 
is [A] ^ g/d, we find that when the production rate, g, is varied while keeping all other 
parameters fixed, r ~ [A]. Otherwise, when the degradation rate, d, is varied while all 
other parameters are fixed, r ~ [A]"^. In general, the switching time is r = T{k,g,d), while 
the population size, [A], of the dominant specie depends on both g and d. Thus, by a 
suitable variation of the rate constants, any desired dependence of r on [A] can be obtained. 
In particular, by increasing k, t can be increased with no effect on [A]. A similar result 
is obtained when g and d are decreased by the same factor. We thus conclude that the 
population size is only one of several factors that affect the switching time. A complete 
description of the switching time should include all the relevant rate constants. 

In Monte Carlo simulations of a switch system with cooperative binding, the switching 
time was found to depend exponentialy on the copy number Q, |l5|. This is consistent 
with the discussion above, but requires a well defined protocol according to which the rate 
constants are varied. 



IV. THE SWITCH WITH BOUND REPRESSOR DEGRADATION 

Consider a different variant of the general switch, in which not only free repressors, 
but also bound repressors are affected by degradation. The bound-repressor degradation 
(BRD) tends to prevent the dead-lock situation in which both A and B repressors are bound 
simultaneously. This is due to the fact that degradation removes the bound repressor from 
the system, unlike unbinding, where the resulting free repressor may quickly bind again. It 
turns out that degradation of bound repressors induces bistability not only at the level of 
the master equation but even at the level of rate equations. 
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A. Rate Equations 



The rate equations that describe the BRD switch take the form 



[A] = g{l - [tb]) - d[A] - ao[A] (1 - M) + 

[B] = g{l - M) - d[B] - ao[B] (1 - [r^]) + a^rs] 
[ta] = ao[A]{l- [r^]) - ai[rA] - ^rK] 

[tb] = ao[B] (1 - [tb]) - ailvB] - rf.M, (19) 

where dr is the degradation rate of the bound repressors. Assuming quasi-steady state for 
the binding-unbinding processes we obtain the Michaehs-Menten equations 

where now k = ao/{ai + dr)- Note that the coefficients of [A] and [B] in the second terms 
in Eq. fl20l) can be considered as effective degradation rate constants. 
For steady state conditins, Eq. (120!) exhibits the symmetric solution 

[{d+drkf + 4:dkg]y^-d-drk 

= = 2dk ■ (2^) 

This solution exists for any choice of the parameters. In addition, in some parameter range, 
two non-symmetric solutions exist. These solutions can be expressed as the solutions of the 
quadratic equation 

ddrk^[A]^ + {gdk + ddrk + dlk'^ - gdrk^A] + gd = 0. (22) 
The condition for the existence of two different solutions of this equation is 

(g - dr) [g{kdr - d)"^ - dr{kdr + d)^] > 0. (23) 

In order for them to be positive the condition g > dr must be satisfied. Thus, the bifurcation 
takes place at 
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(24) 



and the non-symmetric solutions exist for k > kc- Linear stability analysis shows that 
whenever the non-symmetric solutions exist they are stable, while the symmetric solution is 
stable only for k < kc- 

The steady state populations of free A and B repressors vs. k, for the BRD switch, 
are shown in Fig. [61 The results of numerical integration of the rate equations (x) are in 
perfect agreement with the analytical results derived above (solid line). We conclude that 
the degradation of bound repressors induces bistability, even at the level of rate equations. 
The emergence of bistability can be attributed to the fact that the effective degradation 
rate for the minority specie in Eq. (120|) is larger than the effective degradation rate for the 
dominant specie. This tends to enhance the difference between the population sizes and to 
destabilize the symmetric solution for k > kc- 

B. Master Equation 

The master equation for the BRD switch can be obtained from Eq. ([5]) by adding the 
term 



For steady state conditions we find that BRD tends to suppress the peak near the origin 
of P{Na, Nb). For a suitable range of parameters, two separate peaks appear, which qual- 
itatively resemble those obtained for the exclusive switch. However, unlike the exclusive 
switch, there is a narrow corridor with small but non-vanishing probabilities that connects 
the two peaks via the origin. As a result, the switching time for the BRD switch tend to be 
somewhat shorter than for the exclusive switch with the same parameters. The switching 
times, r, vs. the repression strength, k, are shown in Fig. [71 

We now examine in what range of parameters this circuit is indeed a switch according to 
the master equation. Unlike the rate equation where the condition for bistability is clear [Eq. 
(IM]) ]. in the case of the master equation the notion of bistability is more subtle. Thus, in 



dr[Sr^,oPiNA, Nb, Ta + 1, Tb) - Sr^,iP{NA, Nb, Ta, Tb)] + 



dr[5r^flP{NA, Nb, rA, rB + l)- 6r^,iP{NA, Nb, ta, tb)]- 



(25) 
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the analysis below we use the following operational criterion. First we define the two states 
of the switch. The ^-dominated state is defined as the set of all states in which Na > 2 
and Nb — 0, 1. Similarly, the S-dominated state is defined by Nb > 2 and Na = 0, 1. 
The system is considered as a switch if, under steady state conditions, the total probability 
to be in either of these states is larger than 0.99. This leaves a probability of only 0.01 
for all the intermediate states, which the system must visit in order to switch between the 
^-dominated and the S-dominated states. As a result, the switching rate is low. 

We used this criterion in order to find the region in the {k,dr) plane of the parameter 
space in which the BRD circuit exhibits bistabihty. It was found that the BRD switch 
exhibits bistability for large enough values of k, as long as the vahic of dj. is not too different 
from d. U dj./d 1, the process of bound-rcprcssor degradation is negligible and cannot 
eliminate the dead-lock situation. If dr/d ^ 1, proteins bind and quickly degrade. As a 
result, the population of the dominant specie is reduced and bistability is suppressed. 

Within the parameter range in which the system exhibits bistability, we examined the 
dependence of the switching time r of the BRD switch on the parameters g, d, ckq and dr. It 
was found that r exhibits linear dependence on the production rate g and on the repression 
strength k (here, k was varied by changing keeping ai and dr fixed). The dependence of 
T on the degradation rate d was found to be approximately 1/d^. Note that as d was veried, 
we kept dr — d in order that the system remains bistable. Since k depends on dr, it slightly 
varied as well. 

Unlike the exclusive switch, where we managed to obtain the scaling properties of r by 
a simple argument, the BRD switch turns out to be more complicated. This is due to the 
fact that there arc several processes that may lead to the flipping of the switch, such as the 
unbinding or the degradation of the bound repressor. A further complication is that the 
two repressors can be bound simultaneously. As a result, we have not managed to obtain 
an expression for r in the BRD switch. 

V. THE SWITCH WITH PROTEIN-PROTEIN INTERACTION 

Consider a switch circuit which in addition to the mutual repression, exhibits protein- 
protein interactions (PPI), namely an A protein and a B protein may form a complex, AB. 
The AB complex is not active as a transcription factor. 
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A. Rate Equations 



The PPI switch can be described by the following rate equations 



[A] = g{l - [rs]) - d[A] - ao[A] (1 - M) + ai[r^] - ^AB 

[B] = g{l - M) - d[B] - ao[B] (1 - [r^]) + a,[rB] - lAB 
[ta] = ao[A] (1 - M) -aiM 

[r's] = ao[B]{l-[rB])-a,[rB]. (26) 

The parameter 7 is the rate constant for the binding of a pair of A and B proteins. The 
Michaelis-Menten equations take the form 



= Y:;^^-m-'mBi (27) 

For steady state conditions, these equations exhibit a symmetric solution, [A] = [B] , for any 
choice of the parameters. It is the solution of 



-fk[Af + (7 + dk) [A]^ + d[A] -g = 0. (28) 

Since all the coefficients of powers of [A] are positive, this equation has only one positive 
solution. Also, within some range of parameters there exist non-symmetric solutions, given 
by the solutions of 

d-fk[A]^ + {d-f + d^k - g-ik)[A] +6/^ = 0. (29) 

The non-symmetric solutions exist only for the range of parameters in which Eq. ( !29|) 
has two positive solutions. The condition for this can be easily expressed in terms of the 
coefficients in Eq. f l29|) . 

As in the case of the BRD switch, bistability is observed even at the level of rate equa- 
tions. Again, the emergence of bistability can be attributed to the fact that the effective 
degradation rate constant for the minority specie is larger than for the dominant specie, thus 
enhancing the difference between the population sizes [note that the effective degradation 
rate constant for A is (ci + 7[i?]), while for B it is {d + 7[^])]. 
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B. Master Equation 



The master equation for the PPI switch can be obtained from Eq. ([5]) by adding the term 

7[(iV^ + l)(iVB + l)P(iV^ + 1,Nb + 1, r^, rn) - NANsPiNA, N^, ta, r^)]. (30) 

For a suitable range of parameters the steady state solution of the master equation exhibits 
two separate peaks. To draw the range of parameters in which bistability takes place we 
apply the operational criterion used above for the BRD switch. We fix g, d and ai and 
examined the system in the (A;, 7) plane. The results are plotted in Fig. [8] (solid line). 

For small values of 7 (weak PP interaction), the circuit does not exhibit bistability. As 
the interaction strength increases the circuit behaves as a switch for a certain range of 
repression strength k. This range broadens as 7 is increased. Unlike the switch systems 
discussed above, in which the bistability gets stronger as k is increased, the PPI switch is 
bistable for intermediate values of k. This can be understood as follows. Recall that the key 
to the formation of a switch is the elimination of the dead-lock situation. The exclusive and 
the BRD switches deal with this situation directly at the bound repressor level. However 
the PP interaction does not directly affect the bound repressor. To prevent the possibility 
of two proteins bound simultaneously, one of them should unbind and form a complex with 
a protein of the other specie. In order for this to happen, the repressors must not be bound 
too strongly. Therefore, the PPI switch works at intermediate repression strength. As the 
PPI becomes more effective (larger 7) this mechanism applies at larger values of k. 

Enhanced switching properties can be obtained by considering a hybrid system that 
combines PPI and exclusive binding. The resulting switch exhibits bistability in a broader 
range of parameters than the exclusive or PPI switches alone. The master equation for the 
exclusive-PPI switch is obtained from Eq. f|T4|) by adding the term 

-^[{Na + 1){Nb + 1)P{Na + 1,Nb + 1, r^, tb) - NaNbP{Na, Nb, r^, r^)], (31) 

which accounts for the PP interaction. Numerical results, shown in Fig. [8], indicate that 
indeed as expected the exclusive-PPI switch is a better switch than either the PPI or the 
exclusive switch. The parameter range in which it exhibits bistability is broader. Thus, it 
is more robust to variations in the parameters than the exclusive or PPI switches. 
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VI. COOPERATIVE BINDING 

Cooperative binding is found in genetic switch systems such as the phage A switch In 
this case, transcription regulation is obtained only when several copies of the repressor are 
bound simultaneously. This situation can be achieved in two ways. One possibility is that 
repressors bind to each other and form a complex, which then binds to the promoter. The 
other possibility is that the repressors bind separately, but those already bound assist the 
other ones to bind more effectively. In the case of cooperative binding, bistability turns out 
to appear even at the level of rate equations [l7 |. 



A. Rate Equations 

Switch systems with cooperative binding are commonly described by 



l«l = TTW '''' 

where n is the Hill coefficient. It corresponds to the number of copies of the transcription 
factor which are required in order to perform the repression process. Here we focus on the 
case n = 2, and show that these equations exhibit two stable steady state solutions for some 
range of parameters. Imposing [A] = [B] = in Eq. (!32l) . we obtain 



g ~ d[A] - kd[A][B]'^ = 

g - d[B] - kd[B][A]'' = 0. (33) 
Subtracting one of these equations from the other we find that 

- di[A] - [B]) - kd[A][B]m - [A]) = 0. (34) 

Looking for a non-symmetric solution for which [A] ^ [B] , we divide Eq. (IMl) by [A] — [B] . 
We find that /c[A][_B] = 1, or [B] = l/k[A]. Inserting this into Eq. (133!) we get an equation 
for [A]: 
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dk[A]^ - gk[A] + d = 0. 



(35) 



This equation exhibits two distinct stable solutions 



[A] 



gk ± y/g'^k'^ - Ad'^k 
2dk 



(36) 



for k > Ad'^/g'^. This means that the system becomes bistable at the bifurcation point, 
k = Ad'^/g'^. In addition to these solutions, the symmetric solution [A] = [B] exists for any 
choice of the parameters. This symmetric solution is stable for k < Ad'^/g'^ and becomes 
unstable at the bifurcation point. 

B. Monte Carlo Simulations 

Consider a switch system with cooperative binding with n = 2, in which two proteins of 
the same specie bind together to form a complex or dimer. The repression of A synthesis 
is done by dimers composed of two B proteins and vice versa. For example, consider an 
exclusive switch, in which the dimers of A and B cannot be bound simultaneously. To 
account for stochastic effects, we have studied this system using Monte Carlo simulations. 

The rate constant for the formation of dimers is denoted by 7/5. It is assumed that 
dimers cannot dissociate into single proteins, but they can degrade. The degradation rate 
of dimers is denoted hj do- We examined the dependence of the switching time r on all the 
parameters. We found the following properties. The dependence of r on g was found to be 
linear as for the exclusive and BRD switch. The dependence on d is very weak, except for 
the limit in which d is very large. This is because the proteins tend to form dimers before 
they have a chance to degrade. The dependence of r on = ao/«i is found to be well fitted 
by a quadratic polynomial. This means that for sufficiently strong repression, r ~ fc^. 

The dependence of r on the dimerization rate 'jr, [Fig. [9](a)] exhibits interesting behavior. 
For small values of 7d the system is not really a switch, because almost no dimers are formed. 
Therefore the switching time is short. For larger values of 7^) the dimer population increases 
and the system starts to function as a switch. The switching time r increases as the switch 
becomes more stable. But from some point, increasing 7/5 causes r to decrease. This is 
because, very fast dimerization helps the minority specie to form dimers, making it more 
likely to flip the switch. 
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The dependence of r on do [Fig. Mjo)] was found to be well fitted by a cubic polynomial 
in 1/do- This means that in the limit of slowly degrading dimers, r ~ ^/d\). In the limit of 
fast dimer-degradation the system is not bistable, because the population of dimers is too 
small to make the repression effective. 

The switching time for this system was also studied in Ref. where r was presented 
as a function of the average copy number of the dominant specie. However, the copy number 
depends in a non-trivial way on the parameters and cannot be directly controlled. Therefore, 
we believe that in a systematic study of the switching times, it is more practical to examine 
the dependence of r on the parameters themselves. 

Note that there is another important realization of cooperative binding in which the 
promoter consists of two binding sites. When a protein binds to one of them it facilitates 
the binding of another protein to the second site. The effect of this mechanism is qualitatively 
similar to the one shown above for dimers. In general cooperative binding induces bistability 
becuase it forces the minority specie to recruit at least two proteins in order to flip the 
switch. As a result, cooperative binding helps to remove the dead-lock situation in which 
both species are suppressed simultaneously. 



VII. RESPONSE TO EXTERNAL SIGNALS 



Until now our discussion considered only spontaneous transitions between the two states 
of the switch. Here we demonstrate how an external signal may lead to the flipping of the 
switch. In case of the A switch, such an external signal may be, for example, the exposure 
of E. coli infected by phage A to UV light. In the lac circuit, the external signal indicates 
the presence of lactose. We assume that the effect of the external signal is that one of the 
proteins undergoes a conformal change that prevents its binding to the promoter. When the 
signal affects the dominant specie, this may lead to the flipping of the switch. We assume 
that the conformal change is fast and that it lasts for a period of time determined by the 
duration of the external signal. 

We have performed Monte Carlo simulations, where the binding rate ao of the dominant 
specie was set to zero for some period of time (the length of external signal). We calculated 
the probability for a flipping of the switch during 1800 (s), which is roughly the time between 
divisions of E. coli, as a function of the signal length. The results are shown in Fig. [10 
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For short duration of the signal, the switch has a small chance to flip. As the duration 
increases the probability to flip increases too, and so for a long enough signal, the switch 
will eventually flip as expected (the actual switching time depends on the parameters of the 
switch, like the production rate g or the unbinding rate ai. Here we just demonstrated that 
in principle the switch will flip states in response to an external signal). 

VIII. DISCUSSION 

In the rate equations, the meaning of bistability is clear. It typically appears as a result 
of a bifurcation. Below the bifurcation there is a single, stable solution, which becomes 
unstable at the bifurcation point, where two stable solutions emerge. In case of the toggle 
switch, one of these solutions is dominated by A proteins and the other is dominated by B 
proteins. Since both solutions arc stable, the possibility of spontaneous transitions between 
them due to stochastic fluctuations is not included in the rate equation model. 

The objects that participate in regulatory processes in cells, namely genes, mRNAs, 
proteins and promoter sites are discrete objects, and some of them often appear in low copy 
numbers. This, together with the fact that many of the relevant processes such as diffusion, 
degradation as well as binding and unbinding of transcription factors are of stochastic nature, 
requires to consider the role of stochastic fluctuations in these regulatory processes. This 
can be done by using the master equation or Monte Carlo simulations. 

In the master equation, bistability is characterized by two separate peaks in the probabil- 
ity distribution. These peaks should be sufficiently far from each other, with low probabilities 
in the domain between them. As a result, the flow of probability between the two peaks 
is low and the time between spontaneous switching events is long. In order to qualify as a 
switch, the average time between spontaneous switching events must be much longer than 
the time constants of the transcription, translation and degradation processes in the cell. 

For the systems studied here it was found that the general switch without cooperative 
binding does not exhibit bistability bistability either with the rate equations or with the 
master equation. Two other variants, the BRD and the PPI switch systems, were found 
to exhibit bistability both with the rate equations and with the master equation. However, 
the exclusive switch, which is not bistable at the rate equation level, was found to exhibit 
bistability with the master equation. Thus, in case of the exclusive switch it is clear that 
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stochastic fluctuations play a crucial role in making the system bistable. For this system we 
also found an exact phenomenological expression for the switching time in terms of the rate 
constants of the relevant processes. 

Stochastic analysis of genetic networks can be done either by direct integration of the 
master equation or by Monte Carlo simulations. The master equation provides the proba- 
bility distribution of the population sizes of all the mRNA's and proteins in the simulated 
circuit. It can be considered as a distribution over a large number of genetically identical 
cells. The average population sizes and the rates of processes are expressed in terms of 
moments of this distribution. To obtain such distributions from Monte Carlo simulations, 
one needs to repeat the simulations a large number of times and average over them. This 
may be inefficient in terms of computer time, and the statistical errors may be significant. 
On the other hand, unlike the master equation, Monte Carlo simulations enable to follow 
the time evolution of a single cell and directly evaluate quantities such as switching times 
and oscillation periods. 

The number of equations in the master equation set increases exponentially with the 
number of proteins and mRNAs included in the simulated circuit. As a result, the master 
equation becomes infeasible for complex networks. Recently, we have shown that for reaction 
networks described by sparse graphs, one can use suitable approximations and dramatically 



reduce the number of equations 



28|. 



A related circuit, the mixed feedback loop, in which A is a repressor to B and the A and 
B proteins bind to form a complex was recently studied using rate equations |29|, l30| . It was 
found to exhibit bistability within a range of parameters. 



IX. SUMMARY 



Genetic switch systems with mutual repression of two transcription factors, have been 
studied using a combination of deterministic and stochastic methods. These system exhibit 
bistability, namely two stable states such that spontaneous transitions between them are 
rare. Induced transitions take place as a result of an external stimulus. We have studied 
several variants of the genetic switch, which exhibit cooperative binding, exclusive binding, 
protein-protein interactions and degradation of bound repressors. For each variant we ex- 
amined the range of parameters in which bistability takes place. Numerous studies have 
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concluded that cooperative binding is a necessary condition for the emergence of bistabil- 
ity in these systems. We have shown that a suitable combination of network structure and 
stochastic effects gives rise to bistability even without cooperative binding. The average time 
r between spontaneous transitions was evaluated as a function of the biological parameters. 
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FIG. 1: Schematic illustrations of (a) the general switch circuit, that includes two transcription 
factors, A and B, which negatively regulate each other's synthesis; (b) the exclusive switch, in 
which there is an overlap between the promoter sites of A and B proteins, so they cannot be 
bound simultaneously. 

FIG. 2: (Golor online) The probability distribution P{Na,Nb) for the general switch, under 
conditions of (a) weak repression {k = 0.005) where there is one symmetric peak; and (b) strong 
repression (k = 50) where three peaks appear, one dominated by A, the second dominated by B 
and the third in which both species are mutually suppressed. The weights of the three peaks are 
about the same. 

FIG. 3: (Color online) The population sizes of free A and B proteins vs. time obtained from 
a Monte Carlo simulation (a) for the general switch, where the system exhibits fast transitions 
between its three states; (b) for the exclusive switch The bistable behavior is clearly observed, 
where the population size of the dominant specie is between 20 — 60 and the other specie is nearly 
diminished. Failed switching attempts arc clearly seen. The typical switching time is in the order 
of 10^ (s) or roughly 1 day. Bound proteins are also shown. Their fast binding and unbinding 
events cannot be resolved on the time scale that is presented. In both cases, g = 0.2, d = 0.005, 
ao = 0.2 and ai = 0.01 (s~^). 

FIG. 4: (Color online) The probability distribution P{Na,Nb) for the exclusive switch, under 

conditions of (a) weak repression {k = 0.005) where there is one symmetric peak (b) intermediate 
repression {k = 1) where two distinct peaks begin to emerge but are still connected, and (c) strong 
repression {k = 50), where bistability is observed. 

FIG. 5: (Color online) Scaling properties of the switching time r for the exclusive switch vs. the 
protein synthesis rate g, the degradation rate d and the repression strength k. 

FIG. 6: Population sizes of the free A and B proteins vs. k for the BRD switch obtained from 
the rate equations. The parameters are g = 0.05, d = dj- = 0.005, ai = 0.01 and ag is varied. Here 
kc ~ 1.92. Stable solutions are shown by solid lines and unstable solutions by dashed lines. 
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FIG. 7: (Color online) The switching time r vs. k for the exclusive (x), BRD (o) and PPI (A) 
switch systems. The parameters used are g = 0.05, d = dr = 0.005 and 7 = 0.1 (s~^). 



FIG. 8: The range of parameters in the (7, k) plane in which bistability takes place in the PPI 
switch (solid line) and in the exclusive-PPI switch (dashed line), using rate equations (a) and using 
the master equation (b). The other parameters are g = 0.05 and c? = 0.005 (s'^). 



FIG. 9: (Color online) The dependence of the switching time r for the dimers exclusive switch on 
the dimers degradation rate do (a) and the dimerization rate (b). 



FIG. 10: Probability for the exclusive switch to flip during 1800 (s) after the initiation of the 
signal, as a function of the external signal duration. The parameters used were g = 0.2, d = 0.005, 
«! = 0.01 and ag = 0.2 or zero during the signal. 
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